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Abstract 


Approximate  rules  for  evaluating  linear  functionals  are  often 
obtained  by  requiring  that  the  rule  shall  give  exact  value  for  a  certain 
linear  class  of  functions.  The  parameters  of  the  rule  appear  hence  as 
the  solution  of  a  system  of  equations.  Thi3  can  generally  not  be  solved 
exactly  but  only  "numerically".  Sometimes  large  errors  occur  in  the 
parameters  defining  the  rule,  but  the  resultant  error  in  the  computed 
value  of  the  functional  is  small.  In  the  present  paper  we  shall  develop 
efficient  methods  of  computing  a  strict  bound  for  this  error  in  the  case 
when  the  parameters  of  the  rule  are  determined  from  a  linear  system 
of  equations. 
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1.  Introduction 


In  this  paper  we  shall  analyze  mechanical  quadrature  rules  and 
interpolation  formulae  which  have  been  determined  numerically  by  means 
of  solving  a  linear  system  of  equations.  This  process  can  often  not  be 
carried  out  exactly  and  we  want  to  study  the  errors  in  the  computed  value 
of  the  functional  which  hereby  arise. 

In  section  2  we  give  a  general  formulation  of  rules  which  can  be 
found  by  using  the  method  of  undetermined  coefficients  and  outline  a 
computational  process  which  delivers  a  strict  error  bound  in  an  economical 
manner. 

In  the  last  section  we  treat  so-called  Newtonian  feasible  rules 
(See  (7])>  a  class  of  formulae  which  contains  the  Lagrangian  and  Hermitian 
rules  as  special  cases.  These  rules  have  the  pleasant  property  that  they 
can  be  computed  by  a  small  number  of  multiplications  and  divisions.  We 
give  a  general  theoretical  result  on  error  bounds  for  such  rules  and 
illustrate  with  examples  that  it  is  possible  to  solve  problems  in 
integration  and  summation  of  series  in  an  efficient  fashion  by  using  the 
algorithms  in  [71. 


2. 


A  general  class  of  linear  rules 


We  introduce  some  notations  which  will  be  used  in  this  section. 

Let  [a,b]  be  a  closed  bounded  interval  and  let  f;  f^,  fg,  . ..,  l'n 

be  n  +  1  given  functions  on  [a,b]  .  Further,  let  L;  L^,  Lg,  L n 

be  n  +  1  given  linear  functionals  such  that  L(  f ) ,  L^(f),  L^(f'r)  are 

all  defined  for  i  =1,  2,  n  ,  r  =  1,  2,  n  . 

Put  =  L(f^)  r  =  1,  2,  n  snd  let  these  numbers  be  known. 

Sometimes  we  shall  call  yn ,  y_,  ...,  y  moments  with  respect  to  L 
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and  the  system  of  functions  f^,  fg,  ...,  f^  .  This  terminology  is 
motivated  by  the  fact  that  a  wide  class  of  linear  functionals  have  the 
representation 


b 

L(f)  .=  J  f(t)  da(t) 

a 


and  hence 


b 

yr  =  I  rr^  da^  *  r  =  lf  2>  n  * 

a 


We  want  to  approximate  L  by  L  ,  a  linear  combination  of 
Lj,  Lp,  ...,  in  such  a  manner  that  L( f y)  =  L(fr)  ,  r  =  1,  2,  ...,  n  . 
Thun 


(2.1) 


L(f) 


m.L.(f) 


2 


where 


n 

(2.2)  £  m,L.(f  )  =  y  ,  r  =  1,  2,  •  •>>  n  • 

i=l  11  r  r 

We  must  require  that  the  linear  system  (2.2)  has  a  solution.  The  formulation 
(2.1),  (2.2)  applies  for  many  familiar  problems.  We  give  some  examples. 

Example  2:1  A  Lagrangian  integration  rule:  Let  x^,  x . ..,  x^  be 
n  distinct  numbers  and  define  L^(f)  =  f(x^)  ,  i  s=  1,  2,  . ..,  n  and  put. 

b 

L( f)  =  J  f(t)  dt 
a 

Introduce  further  fr(t)  <=  t  ,  r  =  1,  2,  . . n  .  Then 

yr  -  I  t1”1  dt 

a 

Example  2:2  An  Hermit Lan  quadrature  rule:  Let  now  n  be  an  even  number 

and  put  n  «=  2k  .  Select  k  distinct  numbers  x^,  x^,  . ..,  x^  and  put 

Lg^_^(f)  =  f(x^)  >  L^i(  f)  =  f  (x^)  >  i  *  1>  2)  •  ••>  k 

Define  L  ,  f  (and  y^)  as  in  the  preceding  example. 

Example  2:3  A  Lagrangian  derivation  rule:  Let  x  be  a  fixed  number.  Define 

and  f  as  in  example  2:1,  but  put  L(f)  =  f'(x)  .  Then  y^  =  1 

and  yr  =  (r-l)xr-2  ,  r  >  1  . 
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It  is  possible  to  obtain  general  rules  by  appropriate  selection  oi’  the 
function  system  f^,  fg,  •  ••>  f  •  An  obvious  generalization  is  to  replace 
the  interval  [a,b]  by  other  types  of  sets.  Therefore  one  can  extend 
(2.1),  (2.2)  to  the  rules  treated,  e.g.  in  [2]  snd  [9]. 

If  the  coefficient  matrix  of  (2.2)  is  regular  we  can  replace 

(2.1),  (2.2)  with  an  algebraically  equivalent  problem.  Let  namely  A 
be  a  regular  matrix,  n  by  n  and  b,  c,  x  and  u  n-dimensional  column 
vectors.  Then  (2.1),  (2.2)  is  a  special  case  of  the  task:  Evaluate 

m 

y  =  c^x  when  Ax  =  b  .  However,  we  find  immediately  that  we  can  also 
T  T 

write  y  =  b  u  when  A  u  =  c  .  Using  this  observation  we  can  replace 

(2.1),  (2.2)  with  the  alternative  formulation 

(2.3)  L(f)  =  £  cryr 

r=l 


when 


n 

(2.4)  £  crLi(fr)  =Vf)  ’  1  =  2»  n 

r=l 

In  analogy  to  the  usage  in  the  theory  of  linear  programming  we  shall  call 

(2.1),  (2.2)  a  primal  problem,  (2. 5)  and  (2.4)  its  dual. 

We  establish  easily  that  the  duals  of  the  tasks  in  examples  2.1, 

2.2  and  2. 3  consist  of  the  determination  of  certain  interpolating  polynomials. 

We  now  want  to  derive  general  error  bounds  by  using  the  dual  problems 

introduced  above.  Assume  that  (2.2)  has  been  solved  numerically  yielding 

the  approximation  for  rru  ,  i  =  1,  2,  ...,  n  .  Define  Ann  by 

m.  =  m.  +  Am.  and  let  AL  be  the  error  in  L(f)  caused  by  using  m. 
111  —  1 
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instead  of  ,  i  *  1,  2,  . . .#  n  .  Introduce  also  the  residuals 
e  ,  r  -  1,  2,  . n  given  by 

n 

Er  e  yr  "  >  r  s  2#  .  .<#  n  . 


Hence 


(2.5) 


AL 


AmiLi(f) 


when 


(2.6) 


£taiLi<fr> 


®  #  r  B  1#  2>  .««#  n  . 


The  dual  of  this  problem  reads: 


(2.7) 


n 


when 


n 

(2.8)  £  c^L^f^)  =  L^( f)  #  i  ■  1>  2,  *  ••#  n  • 

r=l 

The  formulation  (2.5)#  (2.6)  can  be  used  only  if  the  residuals  are  known 
with  good  relative  accuracy.  This  often  requires  that  they  are  evaluated 
by  means  of  arithmetic  operations  in  a  higher  precision  than  that  which 
was  used  during  the  solution  of  (2.2).  This  drawback  can  be  eliminated 
if  one  uses  (2.7)#  (2.8)  instead.  From  (2.7)  we  get  the  error  bound 
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(2-9) 


|al|  <  e  •  r  where  e  >  m8x|e  |  ,  r  «=  T  |  c  | 

r  r*. L  r 

In  order  to  use  (2.9)  we  need  only  bounds  on  | e r |  and  f  .  The 
latter  quantity  will  later  be  referred  to  as  the  error  factor.  In  the  next 
section  we  will  give  a  theorem  which  expresses  I"  in  terms  of  the  higher 
derivatives  of  f  if  the  rule  defined  by  (2.1),  (2.2)  belongs  to  a 
certain  class. 

Generally  the  error  factor  is  most  easily  found  by  solving  (2.8) 

which  can  be  done  without  too  much  effort.  We  observe  namely  that  the 

coefficient  matrix  of  (2.8)  is  the  transpose  of  that  of  (2.2).  Therefore 

if  we  solve  the  latter  system  by  means  of  Gaussian  elimination  with 

pivoting  we  obtain  a  partition  of  its  coefficient  matrix  which  can  be 

utilized  for  the  subsequent  solution  of  (2.8).  Hence  this  system  can  be 

solved  by  means  of  about  n  operations.  (We  use  the  word  "operation" 

for  a  multiplication  or  a  division.)  Since  the  residuals  can  be  evaluated 
2 

by  means  of  n  operations  the  total  number  of  operations  to  obtain  a 
rule  of  the  type  (2.1),  (2.2)  and  an  error  bound  can  be  written 

(l  +  6/n)(n5/3  +  0(n2)) 


The  error  analysis  can  be  carried  out  in  an  analogous  manner  for  the  case 
when  (2.J),  (2.k)  are  used  instead  of  (2.1),  (2.2). 


3.  Newtonian  feasible  rules 


In  this  section  we  shall  treat  the  case  when  f  is  defined  by 

r 

fr(t)  =  t*”1  ,  r  =  1,  2,  n 

Then  the  solution  of  (2.2)  is  the  coefficients  or  a  polynomial  Q,  of 
degree  less  than  n  .  We  shall  also  require  that  we  can  associate  with 
(2.1),  (2.2)  n  arguments  (not  necessarily  distinct)  in  such  a  manner 
that  Q  can  be  expressed  by  means  of  Newton’s  formula  with  divided 
differences  with  respect  to  these  arguments.  A  rule  meeting  these  conditions 
will  be  termed  a  Newtonian  feasible  rule  (this  definition  is  equivalent 
with  that  in  [7]) . 


Example 

n  =  6  L(f)  =  mjf(O)  +  m?f'(0)  +  rr^T"(o)  +  n^f(l)  +  n^f'U)  +  m6f"(l) 

This  is  a  Newtonian  feasible  rule  since  we  can  introduce  the  six  arguments: 
0,  0,  0,  1,  1,  1  .  If  f  has  two  continuous  derivatives  we  can  express 
these  in  the  form  of  confluent  divided  differences. 

Counter-example 


n  <=  2  L( f)  =  m1f(0)  +  m^f'(l) 
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This  is  not  a  Newtonian  feasible  rule  since  we  need  the  three  arguments 
0,  1,  1  to  express  f(o)  and  f*(l)  in  the  form  of  divided  differences 
but  n  is  only  2.  Still  (2.2)  has  in  this  case  the  unique  solution 
1  ,  =  1  •  We  now  prove  the  general  theorem: 

Let  f  have  n  continuous  derivatives  on  [a,b]  and  let  ( 2 . l) , 

(2.2)  define  a  Newtonian  feasible  rule.  Let  further  the  arguments  associated 
with  the  rule  be  x0,  ...,  x^  .  Define  d^,  >  •  •  •  >  d^  by 

d  =  |f^"'^(t)  |/(r-l)  1  ,  r  =  1,  2,  ...,  n 

1  t€I 

where  I  is  the  smallest  interval  containing  x^,  x0,  ...,  x^  .  If 

c„,  ....  c  is  the  solution  of  (2.1)  then 
1  2  n 

(3.1)  t  |c  1  <  f  d  Y  (1  +  lx  |)  . 

r=l  r  r=l  r  j=l  J 

Proof:  Define  Q  by 

Q(t)  =  ^  crtr_1 

Since  the  rule  is  Newtonian  feasible  we  can  write  Q  under  the  form 

Q(t)  =  r  D  n  (t  -  x  ) 
r-1  r  j=l  J 

where  is  a  divided  difference  with  the  r  arguments  x  ,  x^,  ...,  x^  . 

i  Since  1‘  has  n  continuous  derivatives  there  is  a  number  £  in  I 

sr 

such  that 
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D  =  f  (?  )/(r-l)*  ,  r  =  1)  2)  •  ••>  n 

r  r 


Therefore  the  sum  of  the  absolute  values  of  the  oociT icients  of  y,  is 
less  than  the  sum  of  coefficients  in  Q  defined  by 


n  r-1 


5(t)  =  7*  a  n  (t  +  |x 

rTl  j=l  l] 


But  the  sum  of  coefficients  of  Q  is  Q,(  l)  .  Hence  the  assertion  follows. 
We  observe  that  equality  holds  in  (.5.1)  e.g.  if 

X1  <  0  *  X1  =  X2  =  *•*  =  Xn  and  f(r_l)(-l)/(r-l)’  =  dr  • 

We  conclude  our  analysis  by  discussing  a  few  numerical  examples. 

All  of  these  were  run  on  Stanford's  IBM  %0/bJ.  Its  Algol  W  compiler 

represents  floating  numbers  in  the  form 


,  — .  x 

z  =  x  •  lo 


where  x'  is  allotted  2k  bits  in  single  precision,  5b  bits  in  double. 
Furthermore  x"  is  (if  possible)  so  selected  that  l/lb  <  |x'|  <1  . 

In  all  of  our  examples  we  work  with  Newtonian  feasible  rules.  If 
one  has  to  evaluate  an  expression  in  order  to  get  input  data  such  as 
abscissae  and  moments  this  is  done  in  double  precision.  These  data  are 
afterwards  truncated  to  single  precision.  This  procedure  was  adopted  in 
order  to  insure  that  the  abscissae  and  moments  were  represented  in  full  single 
precision,  independently  of  the  manner  in  which  they  were  obtained. 


The  quadrature  rules  appearing  in  the  examples  were  computed  by  means 
of  the  algorithms  given  in  [7].  The  error  bounds  were  estimated  according 


to  (2.9). 

The  residuals  were  computed  by  means  of  double  precision  arithmetic. 
Thus  they  were  obtained  in  full  relative  precision. 

The  accumulations  to  form  the  scalar  products  which  give  the  computed 
value  of  the  functional  were  done  in  double  precision.  During  this 
computation  the  fact  was  utilized  that  the  product  of  two  single-precision 
numbers  is  delivered  in  double  precision  by  this  particular  machine  and 
compiler. 

It  goes  without  saying  that  a  more  efficient  (but  more  difficult 
to  report)  use  could  have  been  done  of  the  available  resources.  The 
formula 


LL 


li 


|c  E  I 

I  pi 


derived  directly  from  (2.7)  would  presumably  give  smaller  but  still 
strict  error  bounds.  The  computed  value  of  the  error  factor  T 
Indicate  chat  the  total  error  is  bounded  by  a  rather  moderate  multiple 
of  the  largest  residual.  This  could  be  brought  down  most  efficiently  by 
using  double-precision  arithmetic  during  the  evaluation  of  the  weights 
of  the  pertinent  quadrature  rule. 


Example  3:1«  The  integral 


J'  -H 

0  1  +  t 


dt  = 


was  evaluated  by  means  of  Lagrangian  quadrature  rules  with  abscissae 


.1 

« 

(I 


xi  ,  i  *  1,  2,  n  ,  located  in  the  zeros  of  the  function  p  defined 

v  v 

by  g(t)  =  Tn(2t  +  l)  where  Tn  is  the  Cebysev  orthogonal  polynomial 
of  degree  n  .  That  is 


The  integrand  f  is  given  by  f(t)  =  - *  and  hence  the  moments  y 

1  +  t 

yr  B  /  t3>1  dt  =  Vr 
r  0 


In  this  case  the  exact  values  of  the  weights  can  be  computed  by  means 
of  the  formulae  in  [8],  page  127.  We  report  the  following  results. 


Number  of 
moments 

Absolute 
value  of 
observed 
maximum  error 
in  weight 

Absolute 
value  of 
differences 
between 
«/4  and 
computed 
result 

Absolute 
value  of 
largest 
residual 

Error* 

factor 

Estimated 

error 

bound* 

3 

1.2-10"7 

-4 

9.2-10 

1.3-10-8 

1-55 

1.9*10“''' 

6 

3.5-10"6 

4.7-10-6 

2.4- 10'7 

3-24 

7-  •*  10“ r 

4 

Q 

✓ 

1.9-10"3 

2.3-10-7 

1.9*10"6 

5-53 

1.0-1G-' 

1 

■ft 

In  this  and  following  examples  "error"  refers  to  the  error  in  the 
computed  value  of  the  functional  caused  by  the  fact  that  the  weights 
of  the  rule  are  determined  numerically,  not  exactly. 
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The  example  illustrates  the  fact  that  although  the  weights  are  not  very- 
well  determined  the  bound  for  the  contribution  to  the  error  in  the  computed 
value  caused  by  this  may  be  rather  small.  The  circumstance  that  for 
three  moments  the  observed  difference  between  the  computed  integral  and 
n/k  is  larger  than  the  bound  must  be  ascribed  to  the  influence  of  the 
truncation  error. 


1  h+sin  t. 

Example  3:2  Compute  J‘  6  1  1  -  ■  dt 

0  1  +  t 


This  example  illustrates  how  a  suitable  choice  of  a  weight  function  can 
result  in  accurate  quadrature  rules.  These  latter  are  computed  with 
the  algorithms  described  in  [7].  In  this  case  we  take  the  integrand  f 
defined  by 


'f(t)  ■  e 


1 

M-sin  t 


Hence  the  moments  y^,  y^,  yn  are 


yr  =  J 


1  tr~1ln(l/tj 
1  +  t 


They  are  obtained  by  the  recurrence  relation 


*1  *  >  yr  +  yr+l  =  ,  r  •  1,  2,  ...  . 


Lagrangian  rules  with  abscissae  as  in  example  3:1  were  used.  We  give  the 


results 
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i 


I 

r 

r 

r 

r 

r 

r 

r 

r 

r 

r 

i 

i 

i 

i 

i 

i 

i 

i 


Number  of 
moments 

Computed 
value  of 
integral 

Absolute 
value  of 
largest 
residual 

Error 

factor 

Error 

bound 

2 

1.04370 

O 

6.2-10 

1.34 

,,  ,  .  -8 

0 . 3  *  1  1 

3 

1.04362 

5-3-10"8 

1.31' 

7.4.10’8 

4 

1.04362 

6.2*10“° 

1.3  ' 

8.7-10"° 

We  observe  that  the  error  which  can  be  caused  by  inaccurate  weights  is 
neglible  in  comparison  to  the  working  precision. 


Example  3:3 


Evaluate 


■  £  (-D 


r-1 


r=l 


(r‘ 


+  l)^'1 


This  series  belongs  to  the  general  class  of  series  of  the  form 


00 

Z 

r=l 

where  a^  admits  a  representation 

1  r-l 

a  =  J  t  da(t)  ,  r  =  1,  2,  ... 

r  0J  i 

I 

and  the  integrator  a  is  of  bounded  variation  over  [0,1].  a  is  not  i 

dependent  on  r  .  This  fact  can  sometimes  be  verified  by  means  of  a  table 
of  Laplace  transforms  after  making  the  substitution  t  =  e"U  .  Thus  j 

example  3:3  takes  the  form 
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Compute 


da(t) 


when 


J  tr-1  der(t)  = 
0 


We  use  again  Lagrangian  rules  with  abscissae  allocated  as  in  example  1. 
We  report  the  results 


Absolute 
value  of 


Number  of 
moments 

Computed 

sum 

largest 

residual 

Error 

factor 

Error 

bound 

2 

0.47  85  503 

1.2.10"8 

1.41 

1.7.10"8 

4 

0.4819  880 

3.2-10"8 

2.83 

9.1-10-8 

6 

0.4821  010 

5.7-10'8 

4.24 

2.4.10-7 

8 

0.4821  032 

O 

6.010 

5.65 

3.4-10"7 

10 

0.4821  032 

l.J-Kf7 

6.88 

9.5.10'7 

Example  3:4 

00 

Evaluate  s  =  £ 
r=l 

(-1  )r’1e'/'r 

This  series  belongs  to  the  subset  of  the  general  class  discussed  in  the 
preceding  example  where  the  numbers  a^  introduced  there  are  the 
moments  of  an  nondecreasing  integration  a  .  That  is  we  can  write 

s  ■  /  Sft  ao,(t) 
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when 


J*  t'"1  dar(t)  =  eVr  ,  r  =  1,  2,  ... 

0 

and  in  addition  at  . 

As  shown  in  [l]  the  integral  is  bounded  between  the  values  which  result 
if  certain  quadrature  rules  of  Gaussian  type  are  applied  to  the  integral 
for  s  .  These  rules  can  be  computed  by  algorithms  given  in  [5],  [3], 
[4],  and  [6],  One  must,  however,  solve  a  non-linear  system  of  equations 
which  is  unique  for  each  series.  This  is  avoided  in  the  following  way: 
Since  the  derivatives  of  the  function  f  defined  by 

f(t)  =  (1+t)-1 

do  not  change  sign  we  can  by  means  of  Kermite-interpolation  construct 
two  polynomials  P..  and  P_  such  that 

X  d 

?1(t)  <  f(t)  <  P2(t) 

Let  ti  be  the  number  of  moments  needed  in  the  quadrature  rules  and  let 
xi  and  ti  be  distinct  points  in  (0,l).  The  conditions  defining  P^ 
and  P0  must  be  of  the  form: 

n  =  2k+l 

PjU)  =  P2(o)  -  f(o) 

P^t.)  =  f(tA)  P2(^)  =  f(x.)  i  -  1,  2,  ...,  k 

pjU^  =  f,(ti)  p^x^  =  r’(x±) 
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n  *  2k 


P^)  =  f(tt) 
p^tj)  =  r,(ti) 


P2(xi)  =  f(xi) 
PJ(xi)  =  f'^) 
P2(o)  =  f(0) 
P2(l)  =  f(l) 


1  —  3.  y  III)  lc 


Thus  we  only  have  to  solve  linear  systems  of  equations  (which  can  be  done 

o 

with  n  operations).  Since  and  do  not  depend  on  the  moment 

sequence  the  same  P^  and  P^  applies  to  all  series  of  the  class 

discussed  here.  In  example  3:3  we  computed  upper  and  lower  bounds  by 

v  v 

allocating  the  interior  touching  points  in  the  zeros  of  Cebysev 

v  v 

T^(2x+1)  where  is  a  Cebysev  polynomial  with  degree  equal  to  the 

number  of  such  points.  We  give  the  results 


Number  of 
moments 

Difference 

between 

bounds 

Estimated 
error  in 
upper 
bound 

Estimated 
error  in 
lower 
bound 

Difference 

between 

computed  bounds 
Gaussian  rules* 

3 

5.7-lCf5 

1.3-10"7 

1.1* 10"7 

2.5-10"5 

h.9-10”5 

3.  o«lo”7 

3.2-10"7 

9.9-10"b 

9 

3.0- lO-7 

1.7-10"'" 

6.9-10"7 

U.o-icf8 

For  this  computation  double  precision  was  used  throughout. 
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